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Abstract 

The aim of this paper is to investigate new numerical methods to compute travelling 
wave solutions and new ways to estimate characteristic properties such as wave speed for 
stochastically forced partial differential equations. As a particular example we consider the 
Nagumo equation with multiplicative noise which we mainly consider in the Stratonovich 
sense. A standard approach to determine the position and hence speed of a wave is to compute 
the evolution of a level set. We compare this approach against an alternative where the wave 
position is found by minimizing the norm against a fixed profile. This approach can 
also be used to stop (or freeze) the wave and obtain a stochastic partial differential algebraic 
equation that we then discretize and solve. Although attractive as it leads to a smaller domain 
size it can be numerically unstable due to large convection terms. We compare numerically 
the different approaches for estimating the wave speed. Minimization against a fixed profile 
works well provided the support of the reference function is not too narrow. We then use 
these techniques to investigate the effect of both Ito and Stratonovich noise on the Nagumo 
equation as correlation length and noise intensity increases. 



1 Introduction 

The effects of stochastic forcing on solutions of stochastic partial differential equations (SPDEs) 
have recently received a great deal of attention in appUcations ranging from material science, 
atmosphere modelling to neural science. Of particular interest are the effects of noise on travelling 
waves and fronts as these are often physically important solutions. We use the term travelling wave 
to include fronts and waves and develop in this paper numerical methods to solve for stochastic 
versions of these objects. 

We consider the stochastic PDE 

dt + g[u,t)odW{t), given u(0) = a; G K. (1) 

Although the noise term is written here as a Stratonovich integral we consider below the noise in 
both a Stratonovich and Ito sense. For ease of exposition we take u : M x M+ — >■ M although similar 
numerical proceedures have been applied to systems of PDEs and waves in K^, see for example [5]. 
For additive noise the function g is taken as a constant whereas for multiplicative noise g depends 
on u. In the case of no noise [g = 0) we recover the deterministic PDE 

Ut = Uxx + f{u), given u(0) — u", a; € M. (2) 

In the deterministic case the analysis of travelling waves both analytically and by numerics is a 
mature field. This is not the case for SPDEs where much of the analysis is performed for specific 
equations or for the case of small noise. Indeed with stochastic forcing existence for all time of 
these waves is not guaranteed and the definition of quantities such as wave speed vary from system 
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to system. Typically the position of a stochastic travelling wave is determined from the position of 
a level set. For small noise the centres of these fronts can be shown to follow a rescaled Brownian 
motion, see [331 |5J [Uj . In the case of multiplicative noise the front may exist for all times and the 
wave front may have compact support |36j and is well-defined over some time varying interval in 
space [a{t),b{t)] and takes stationary values out side of this interval. There is a well developed 
literature for the stochastic Fisher-Kolmogorov-Piscounov equation [12l [TTJ [26l [39] with waves 
defined in this form. Multiplicative noise is seen to change the wave speed of the wave and the 
position is seen to diffuse from the mean (or Goldstein mode) , for reviews see [121 [27] ■ However, 
it is not our aim to replicate these results here. We investigate different ways to measure the 
wave speed, using the level set approach as well as a new approach of minimizing the norm 
of the wave against a fixed profile. Furthermore we apply different computational techniques to 
compute time dependent waves - this inlcudes freezing the wave to stop it from travelling. We 
extend a numerical method introduced in [6 for deterministic PDEs of the form This method 
freezes the wave in the computational domain by adding a convection term to the equation to 
compensate for the movement of the wave. The convection term that gives the speed of the wave 
is determined from an extra algebraic condition from the minimization and the wave speed 
is explicitly solved for as a time dependent quantity. Convergence of this method for PDEs was 
considered in [37] and stability of the wave considered in [55]. We use these techniques to obtain 
new computational results on the effect of multiplicative noise on the wave speed of the Nagumo 
equation with Stratonovich and Ito noise. In contrast to |161 [31 [5], we present new numerical 
results on the effect of spatial correlation length on the width of the wave and compare different 
measures of the wave speed. These results incude our new idea of comparing the wave to a 
reference function. Furthermore we present new results for Ito and additive noise. 

Numerically we solve for the wave profile and a time dependent wave speed for ([!]) which 
is, in the case of stochastic forcing, a random variable. As a specific example to illustrate the 
computational method and to compare against existing techniques we consider the scalar Nagumo 
equation [18] 

du = [uxx + u{l — u){u — a)] dt + [v + /iii(l — u)) o dW. (3) 

With 1/ 7^ and ^ = we have additive noise and /i 7^ the noise is multiplicative. For 
multiplicative noise we have that u = and u = 1 are stationary and numerical simulations 
suggest a wave exists between them [S] [2] ■ The deterministic equation 

ut = Uxx + ""(1 - u){u- a), u{x, t) e M, a; e M, t > 0, (4) 

is often used for testing algorithms since travelling wave solutions u{x, t) — Udct{x — ct) connecting 
the stationary points u_ — 0, U-^- = 1 of this equation are explicitly known besides other explicit 
solutions, such as pulses, sources and sinks Tl[9]. These travelling wave solutions depend on the 
nonlinearity and the leading profile of initial data u". Define the function Uk{x) by 

ukix) = {l+e-''^y\ (5) 

We use this function to specify both initial data w° and reference functions that have different 
profiles (by varying k). For a G (0, 1/2] there is a unique asymptotic travelling wave where as for 
a £ (—1, 0] the asymptotic profile and the wave speed depends on the leading profile e''°^ of the 
initial data as a; — >■ (X). We summarize results below for the deterministic Nagumo equation, these 
are found, for example, in |16j . 

• For a e (0,1/2] the solution u = Uk, with k — l/\/2 is asymptotically stable and all 
initial front data is attracted to this wave. The asymptotic wave speed is given by 
c=-V2 (i - a). 

• For a £ (—1/2, 0] if the initial data u° — Uk,-, has ko > k^, — —aV2 then the asymptotic speed 
is given by c = — \/2 (^ — a). If ko < fc* then the asymptotic wave speed is > (fcg — a)/kQ. 

• For a € (—1,-1/2] if the initial data = Ukg is such that ko > k-^ = then the 
asymptotic speed is given by 2k^ . If fco < k^ then the asymptotic speed is > 2k^ . 
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The outHne of the rest of the paper is as foUows. We review in Section [2] the computation of 
travelhng waves in the deterministic case and extend to the case of stochastic forcing. We discuss 
measures of wave speeds of a stochastic travelhng wave and discuss the numerical approximation. 
In Section |3] we illustrate the numerical method on the Nagumo equation with multiplicative 
noise. We compare solving the stochastic partial differential equation (SPDE) and the stochastic 
partial differential algebraic equation (SPDAE) where the wave is frozen by minimizing the 
distance between a reference function and the travelling wave. We compare the different measures 
of wave speeds. We illustrate from numerics that although a feasible method it can lead to 
numerical instability. We investigate the effect of the choice of reference function u in Section [3.2.2[ 
In Section |3.3| we present new numerical results for Ito and Stratonovich multiplicative noise on 
how the wave speed changes with noise intensity. We alos present new results on the effect of 
the spatial correlation length. Additive noise is considered in Section [3.5| where we again examine 
the wave speed with noise intensity and illustrate the SPDAE approach when new travelling 
waves are nucleated. Finally we consider weaker versions of the stochastic travelling wave fixed in 
the computational domain by mean wave speeds and then discuss the results and computational 
method. 



2 Stochastic travelling waves 

In this section we introduce the (stochastic) differential algebraic equations that we use to define 
the travelling wave problem. We start by reviewing the more familiar deterministic case before 
considering the case with stochastic forcing. In both cases we reduce the infinite problem to finite 
dimensions by truncating the computational domain and discretizing in space. 

2.1 Deterministic PDE and discretization 

Let us assume that equation ^ has a travelling wave solution u, so that u can be written as 

U{x,t) = UdctiO^ ^^X~Xdctt, (6) 

where Udct G C^(K,M™) denotes the waveform and Adct its wave speed. In a comoving frame 
v{(,t) — — Xdctt,t) equation ^ reads 

Vt=v^^ + Xd,tv^ + f{v), ^eR, t>0 (7) 

of which the travelling wave Udet is a stationary solution. Since the wave speed Adet is generally 
unknown we transform equation ^ into a co-moving frame with unknown position j{t), i.e. we 
insert the ansatz v{x,t) — u{x — ^{t),t) into Then we obtain 

Vt = Vxx + + f{v), (8) 

where A(t) = j'it). In order to compensate for the additional variable A we add a so called phase 
condition 

= i^{v, A) (9) 

which together with ([s]) forms a partial differential algebraic equation (PDAE) [B]. The position 
7 of the wave can then be calculated by integrating 7' = A, 7(0) = to get 

7(t) = / X{s)ds. (10) 
Jo 

For the numerical implementation we need to truncate the spatial domain from x €M.to x d [0, L] 
and impose appropriate boundary conditions such as Neumann, Dirichlet or projection boundary 
conditions |37l. We then solve (|8| and ([9| for x S [0, L]. In contrast to traditional determinis- 
tic travelling wave computations where the steady states of Q are solved for with appropriate 
boundary conditions this method does not rely on A being a constant wave speed. 
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Thus far we have not discussed the choice of the phase fixing function -0 in (|9]). Since the 
phase condition only selects one representative out of the infinite family of solutions, there is some 
freedom of choice here. The simplest phase condition is to align the solution with respect to a 
given reference function u. It is natural to take the norm since we then the minimum can be 
found by differentiation and equating to zero (see eg [7]). For two functions v and w to minimize 
over shifts in space y the norm: miuj^ W'^'i^j ') ~ — y, ■)\\2- Differentiating and equating to 
zero we find that 

(t;(x, •) — iu{x — y, ■))wxdx = 0. 



So for the PDE we take the phase condition to be : 

V'fix(w) = {Ux,U-u). 

This choice was termed the template fitting method in [32J. 

In our numerical simulations we will discretize in space using standard uniformly spaced finite 
differences, so we discretize on a finite grid xq, Xm, u{xj) = Uj. For the second derivative with M 
points and spatial step Aa; we approximate the derivative dxx ~ A where A = -^^B G ]^m-2,m-2 
for Dirichlet boundary conditions and for Neumann boundary conditions, 



-2 2 
B 
2 



A = ^ I B I e M*^>^ with B = 



(-2 1 \ 



'■■ '■■ 1 

V 1 -2 / 



We choose not to use periodic boundary conditions since we compute a travelling front rather than 
a pulse and this would require a domain of twice the size and also introduces two travelling waves 
that travel in opposite directions. Neumann and Dirichlet boundary conditions were shown to work 
well in the deterministic case in [371 EH] ■ For the first spatial derivative we introduce {Diiu)j = 
{uj+i - u.j)/Ax, {DLu)j = {uj - Uj-i)/Ax, {Dcu).j = (uj+i - iij_i)/(2Aa;) for j = 1, ...,Af - 1 
using Dirichlet boundary conditions — um — 1r or Neumann boundary conditions — 
Ui,um-i = w^v/. For convection terms we either use Dc or, where up- winding is an issue, we 
choose the appropriate D^, or a weighted combination [B] 

9, « D,, g-^^Dl + (1 - e'P'')DR, (11) 

where /3 is a parameter (/3 = or /3 = ^ in our computations) and in what follows h will be some 
function of the wave speed. Recent work by Hairer and Voss [17] examine the discretization of the 
advection term uu^ for the stochastic Burger's equation and show that the limit is dependent on 
the discretization. The form of advection for the Burger's equation and considered here is different 
and as we can compare to cases without the advection term we did not note and discretizatoin 
dependent 

Discretizing in space with N grid points and after eliminating the boundary conditions we 
obtain the following DAE system for A e M and v e M^^^ for Dirichlet or u e for Neumann 
boundary conditions 

v' ^ Av + X{Dxv + Tj)+ fiv) + ^ ^^^^ 
= (ux, V -u), 

where the vectors ip, rj are used to deal with the boundary conditions. This system can be solved 
by using appropriate DAE solvers [1] or we can use a linear implicit Euler method to obtain the 
fully discrete system 



(13) 



= {Dcu,v''+^ 



u) 
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which leads to 



AtA 



Ax Dc-''^ 



-At{Dx^v''- 




AtU{v^) + ^) 
{Dcu, u) 



Note that for the reference or template we use the central difference approximation Dc since 
this is most accurate and convection instabilities are not an issue for this term. 

Under a uniqueness assumption of the travelling wave of the PDE it was shown in [37 that for 

L- 



oo and Ax — ^ the stationary solution of ( 12 ) converges to the exact travelling wave solution. 
Moreover the solution of (12 1 inherits the nonlinear stability properties of Numerically we 
observe below that the DAE system correctly computes the travelling wave depending even when 
the travelling wave is not unique, see Section [3] 



2.2 Stochastic PDE and stochastic travelHng wave 

We seek travelling wave solutions to the Stratonovich SPDE 

du = [uxx + f{u)] dt + g{u) o dW, u(0) = u' 



(14) 



or the ItoSPDE 

du = [uxx + f{u)] dt + g{u)dW, u{0) = u° (15) 

with g(u) = V + ^ih{u), where v and fi are parameters that allow us to consider additive and 
multiplicative noise. Results on the existence of a solution for (14) and (15) with x G K domain 



can be found in [jOj. For the stochastic Allen-Cahn equation with spatially smooth, bounded 
additive noise existence is shown in [21] . A recent paper [33] shows existence for the non-Lipschitz 
cases for space time white noise. 



We truncate the infinite domain and consider (14) (or (15)) on a large finite domain so that 
X £ [0, L] with either Neumann or Dirichlet boundary conditions. For the finite domain with 
X e [0, L] we refer to [TD] Theorem 7.4] with / and g satisfying global Lipschitz conditions, [55] for 
some weaker conditions and recent results in |19j for non-Lipschitz g. For the stochastic Nagumo 
equation we have f{u) = u{l — u){u — a) and take h{u) = u{\ — u). 

We consider noise W{t) to be a Q- Wiener process [lOj, and assume that covariance operator 
Q and the linear operator dxx have the same eigenfunctions If the covariance operator has 
eigenvalues C,j > then we can write 



(16) 



for independent Brownian motions /3j. We take space-time noise that is white in time with 
exponential decay in the spatial correlation length ^ > in which case 



E{dW{x,t)dW{y,s)) = C{x- y)S{t,s), C{x) = —exp^ 



1 

2^ 



TTX 

'4^ 



We approximate using (16) by taking C„ = exp(— S-^-i), where \j = and L is the length of 



the interval 135, "^T 



In the Stratonovich case ( 14 1 we can eliminate the systematic effects of the noise on the drift 



and convert to an Ito integral. This gives an additional term to the nonlinearity so that the 



Stratonovich SPDE (14) is equivalent to the Ito SPDE 



du 



+ f{u) dt + g{u)dW 



(17) 



where f{u) = f{u) — C{0)g' {u)g{u), see for example [THj. We can also convert from the Ito 



interpretation ( 15 ) to a Stratonovich by 



du 



Uxx + f{u) dt + g{u) o dW 



(18) 
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where now /(u) — f{u) + C{0)g'{u)g{u). 

We discretize the SPDE in space by finite differences and evaluate the noise on the spatial grid, 
fn time we discretize with a constant time step At. For the noise term we have an increment 

3 el- J.J] 



where £,j ^ N{0, At). To compute directly with the Stratonovich noise for ( 14 ) we use the standard 
Heun method [151 HO] and also the semi-implicit Euler-Hcun method 

z = u" + .g(u")Aiy„ 

=u" + At [Au"+' + /(^.") + (giz) + giu''))AW^ ^^^^ 

where AWn in an increment of the noise and if arises from the boundary conditions. 

Although intuitively it is understood what is meant by a stochastic travelling wave it is not 
easy to find a definition in the literature, however for a review see [121117]. Typically a stochastic 
travelling wave and speed is either defined by the evolution of a level set such as in [IHl [33] or 
through the evolution relative to a deterministic wave, such as through a small noise expansion 
such as in 23 . We will apply our methods to the case where in the deterministic case the travelling 
wave is known to be unique and also where it is not unique. 

Consider the SPDE with a well defined wave with compact support as defined in [36 , so that 
at u{—oo, • ) = u(oo, • ) — We can then define a travelling wave and wave speed using the 
points 

a{t) :— supjz : u{x,t) = < z}, b{t) :— sup{z : u{x,t) = u+,x > z}, (20) 

and in addition we can take the 'mid point' level set of a wave 

c(t) := sup{z : u{x, t) — (m_ + u^)/2, x < z}. (21) 

These level sets define the position of the travelling wave. Note that we take the supremum as 
there may be multiple crossings through the level set (see for example Figure [s]) for the front. 
Given the positions an 'instantaneous' wave speed can be determined from the a, b and c by 
differentiation. We report below a wave speed Az{t), z e {a, b, c} defined by 

A,Xt)=E^^^^^ ze{a,6,c} (22) 

where the expectation is taken over the number of realizations. We may choose the initial time to 
to either be at the start of the computation tQ — or some later time (ip > 0) to avoid transient 
effects. We believe This differs from the definition of wave speed used in computations by [21 HI Hi] 
where they report z{t)/t for z G {a, 6, c}. Numerically, these level set points a,b,c are found by 



evolving the SPDE ( 14 ) directly and interpolating over the grid. 

If we assume the wave has some long time invariant speed an alternative definition of the wave 
speed is to fit a linear polynomial -PAfit to the data (t,¥,z{t)), z € a,b,c, t > tQ where 

Pa,,M ■■= ^fitt + K (23) 

Wave speeds may then be estimated by Agtzj ^ e a,b,c where we may take <o > to avoid 



transients. Although this is a trivial extension of wave speed defined by ( 22 ) we have not seen it 
reported in the literature. 

Finally we introduce a novel measure of the wave speed fpr SPDEs through minimization of 
the norm ||u(a:,t) — u{x — y,t)\\'^ against a fixed profile u. This is similar to the freezing the 
wave in the deterministic case. We solve the SPDE and compute the position j{t) of the wave. 
We then move the reference function relative to the travelling wave solution u. That is we solve 
SPDE 



du = [uxx + fiu)] dt + g{u) o dW, u(0) 



,0 



u 
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and couple this to a reference function u{x,t) that moves so that 

{u^{t),u{t)-u{t)) =0. (24) 

NumericaUy this requires interpolation onto the spatial grid at each time step. We compute an 
instantaneous wave speed X{t) and this is related to the position of the wave through j{t) = 
Jq X{s)ds. A wave speed Amin is then defined through the time average of the instantaneous wave 
speed A 

A„,in(t) = / E{\{s))ds. (25) 

t - *0 J to 

So far we have not commented on the choice of profile ii to minimize against. 

If a deterministic travelling wave profile exists then this is a natural choice for ii. In examples 
where we do not have an analytic expression for the deterministic travelling wave ii then this can 
be solved for simultaneously or a sample solution solved for and saved. However, this is a matter 
of choice and we could minimize the norm against any fixed profile. One obvious choice of 
profile ii is to take the initial data, so that u — u{0). Note that the choice of u is important. In 
particular we illustrate in Section |3.2.2| that for a given u the minimization may not be unique 
and may fail numerically if u has small support. 

We again note that from the position data {t,E,^(t)) we can fit a linear polynomial 

PA^it) := A^t + K (26) 
to obtain an alternative estimate of the average wave speed. 



2.2.1 Freezing the stochastic travelling wave 



Inspired by the deterministic fixing of a wave we freeze a stochastic travelling wave relative to a 
reference function u, this allows direct computation of the wave speed from the minimization. 
We take u to be a fixed continuous function and we evaluate it numerically at the grid pouints. 
We examine the Stratonovich SPDE in a co-moving frame as we did for the deterministic case. 
First let us examine the effects of a shift in space on the covariance of the noise noting that 

E{dW{x + r{t),t),dW{y + r(s), s)) = C{x + r{t) - y - r{s))5{t ~ s). 

We see that for noise that is white in time the covariance is the same in the two frames. 

Let's consider the Stratonovich SPDE transforming to the moving frame m(x, t) — v{x + ^{t), t) 
we have 

dv{x + 7(i), t) — dj{t) + dv 

and so formally we can write 



dv 



, d^jt) 

^XX T "^X 

at 



fiv) 



dt + g{v)odW, v(0) 



We introduce the phase condition to determine c?7 and minimize the norm, i.e. {ux, v — u) ~ 
If we define the random variable A so that 7(t) ~ \{s)ds then we seek to solve the SPDAE 

dv = + \v^ + f{v)] dt + g{v) o dW, v{0) = u° 



[Vxx 

= (Ux,V - u) 



(27) 



In order to compute directly with the Stratonovich noise for ( 14 1 we use either the standard Heun 
method [HI [50] or the semi-implicit Euler-Heun method 



z 

,n+l 



giu^)AWn 

At [Au''+^ - 



A" (i?A"K" + ry) +/(?.")+ ^] +-(5(z)+5(u"))AVF„ (28) 
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The algorithm yields an approximation m", n ~ 0,1,2, .. . to u{nAt) and A", and approximation 
to A(s). This gives us a numerical scheme for the SPDAE in which the stochastic travelling wave 
is frozen. 

We have a time-dependent random variable A(t) that we call the instantaneous wave speed. 
Of more physical interest is the time average of this quantity A^^^ that we call the wave speed 
and report 

^{Z-Ar t^Hs)ds. (29) 



mm J. J. 

t - to J to 

The instantaneous wave speed A gives the position 7(f) of the wave. If we assume the wave has 
some long time invariant speed this can be estimated from j{t) by fitting a linear polynomial 

Pj^f.4t):=Ai^'^t + K (30) 

to the data from freezing the wave {t,E'y{t)), for t > to >0. 

The literature on solving stochastic DAEs is in its infancy, however there are some analytic and 
computational results mainly arising from examining noise in circuit simulations, see for example 
[331[3nilllllll]- We are not aware of work on the existence directly for SPDAE. 



Computing a travelling wave through (27) or we introduce the random variable A which is 
used to freeze the wave. We could however define a weaker versions by taking statistics of A. For 
example we can take the time-averaged wave speed A(i) for each realization 

dv=^[v^^+A^i^v^ + f{v)]dt + g{v)odW, v{0) = u° 
= V(w,A). ^ ' 

Other weaker forms of travelling wave solution are possible where the instantaneous wave speed 
A or time average wave speed A of an individual realization is replaced by its expectation over 
realizations, for example 

dv^[v.,,+EiX)v, + fiv)]dt + giv)odW, viO) ^ u° 
= tP{v,X); 

and 

dv = + E (A) v., + fiv)] dt + g{v) o dW, w(0) = u° 
= ^(«,A). 

Using the sample mean of A and A for fixing we are essentially using a "group velocity" to fix 
the wave and as a result the mean profile will contain a spread as each individual realization is 
not fixed at the same point. By taking these weaker notions of wave speed to freeze the wave we 
observe spreading of the front profiles, as discussed in P^ . 



3 Results for the Nagumo Equation 

We compare the different estimates of the wave speed and apply the technique of freezing the 
wave to the Nagumo equation ([s]) for both multiplicative and additive space-time white noise. For 
the majority of our simulations we take Aa; — 0.1, At — 0.05 and a spatial domain of i = 500 
or L = 800 with Neumann boundary conditions and integrate till t — 100. We compute 100 
realizations simultaneously. In our computations of wave speeds unless stated we take to = to 
reduce transient effects and we drop the dependence of the computed wave speeds on the time t 
and so report A^in, A;^'^^, A^, A^ and A^t^, z e {a, b, c}. 

3.1 Deterministic PDE 

Before we examine the stochastic PDE we briefly examine deterministic computations. We point 
out some features of computing the travelling wave and speed by direct simulation of the PDE 
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k = 1/v^ 


Theory 


■^min 


Aa Afc Ac 


Afit a Afitf, Afitc 


fco = 1/^2 
fco = 0.1 


1.06066 
> 2.6 


1.06047 
2.59741 


1.06025 1.06027 1.06026 
2.59690 2.59689 2.59689 


1.06025 1.06026 1.06026 
2.59689 2.59689 2.59689 



Table 1: Different measures of the wave speed computed from solving the deterministic PDE. To 
compute Amin the profile u travels to minimize the LP' norm ( 24 1 
are from ( 22 ) the level set and Agt 



the wave speed agree to 4 decimal places. 



from ( 23 ) from fitting, 



Estimates of wave speeds A^, 
z e a, 6, c . We see these measures of 



k -- 


= l/V2 


Theory 


mill 


Aa 


Afc 


Ac 


Afit a 


Afit 6 


Afitc 


ko 


= 1/V2 


1.06066 


1.06052 


2.0e-07 


-l.le-06 


1.5e-08 


1.8e-07 


-9.5e-07 


1.4e-08 


fco 


= 0.1 


> 2.6 


2.60048 


3.1e-06 


-3.6e-07 


-1.4e-07 


1.3e-06 


-8.4e-08 


-6.0e-08 



Table 2: Wave speeds computed from solving the PDAE and freezing the travelling wave. The 
fact that the wave does not move in the domain can be seen from the level set wave speeds A^ 
and Afit^, z G {a,b,c} which are close to zero. The wave speed A^'^^ agrees with that computed 
for the PDE given in Table [ij 



([2]) versus freezing and solving the PDAE (12). In particular we examine the regime where the 
travelling wave is not unique and the theory of |37j on freezing the deterministic case no longer 
holds. For a = —0.25 the asymptotic travelling wave and wave speed depends on the leading 
profile data of the initial data We take two initial profiles u"(a;) = Uko{x) with fco — 0.1 < fc* 
and fco = l/\/2 > fc*. To compute the speed by minimization we present results with reference 
functions u — with fc = ^/2. Results with a reference function with fc = 0.1 are identical (see 
also Section 3.2.2 for comments on the choice of reference function). 

We show in Table [l] wave speeds computed from direct simulation of the PDE. For the PDE 
Amin is computed by moving the profile u at the computed wave speed from the minimization 
using the condition ( [24| . There is good agreement between the computed wave speeds, although 
Amin appears to have converged faster than the other measures of the wave speed to theoretical 
value. When we freeze the wave in the computational domain and solve the PDAE we see from 
Table [2] that the wave is frozen (to single precision) since the level set positions of a{t),b{t) and 
c(t) do not change and hence wave speeds Agt^, z G a,b,c from fitting the linear polynomial 
are zero (to single precision). The wave speed A„ii„ estimated from freezing the wave and the 
minimisation the norm agrees with the wave speeds computed from the PDE (and is in fact a 
better approximation to the theoretical values). 



3.2 Stochastic travelling wave and frozen wave 

To illustrate computations for the stochastic PDE we start by taking Stratonovich multiplicative 
noise with /i = 0.1 and a correlation length of ^ = 0.1. In Table |3] we show results from solving 
the SPDE with the same single realization of the noise using the same two different sets of initial 
data and two difference reference functions as for the deterministic case. Since we have taken 
the same noise realization, when initial data is the same our measures of the wave speed A^, and 
Afitz, z £ a,b,c are identical and independent of the reference function u. The choice of reference 
function does change the wave speed measured by the minimization in approximately the fourth 
decimal place (compare Amin or A.^ for the different fc values in Table [3|. This small difference is 
due to a combination of interpolation errors and is not seen for the SPDAE below where we do 
not need this interpolation. If we compare the values of the deterministic PDE Table [T] and SPDE 
case Table [3] for the single realization we see the wave speeds with noise are slightly larger than the 
deterministic case. In Figure ll] we plot the result of a single realization in (a) for the SPDE with 
initial data vP with fco = l/v2! The wave front is initially a,t x ~ 200 and travels to x ~ 500. For 
the two different initial data we have plotted in (b) the two distributions of the instantaneous wave 
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to 


= 50, ti = 


100 








Ab 


Ac 


Afita 


Afitb 


Afitc 


ko 


= 1/V2,k 


= 1/V2 


1.07575 


1.06965 


1.07484 


1.06746 


1.07536 


1.07149 


1.07147 


1.06935 


ko 


= 1/V2,k 


= 0.1 


1.07538 


1.06987 


1.07484 


1.06746 


1.07536 


1.07149 


1.07147 


1.06935 


ko 


= 0.1, fc = 


0.1 


2.77250 


2.78015 


2.80594 


2.74516 


2.76020 


2.79657 


2.72070 


2.78112 


ko 


= 0.1, fc = 


l/%/2 


2.77146 


2.78178 


2.80594 


2.74516 


2.76020 


2.79657 


2.72070 


2.78112 



Table 3: Wave speeds computed from solving a single realization of the SPDE with noise intensity 
fi = 0.1 and correlation length ^ = 0.1. To compute Ay^in the profile u travels with the appropriate 
speed found by minimization of the norm. 



speed A used to compute the wave speed through the minimization (with k = l/\/2). The mean 
of these distributions gives the corresponding wave speeds, 1.07575 for ko = l/%/2 and 2.77146 for 
fco = 0.1. For initial data fco = l/\/2 with a wave speed of 1.07575 the variance of A is smaller. 
In (c) we have plotted for the two different initial data sets the instantaneous wave speed X{t) 
and the corresponding time averaged wave speeds Amin(i) with to = 0, ti — t. We see faster 
convergence of the wave speed for initial data fco = l/%/2 and again the reduced variability in the 
instantaneous wave speed. 



(a) (b) (c) 




Figure 1: (a) Space-time plot of a single realization of the SPDE showing a travelling wave. In 
(b) distributions of the instantaneous wave speeds A computed for the two different initial data 
sets. In (c) we plot A(i) and the time averages A,„in(i) with <o = 0, ti = t. 



Let us now compare to a single realization where the stochastic travelling wave is frozen and we 
solve the SPDAE (27 1 using the Heun method (28). Results on the wave speeds for the two initial 
data sets and two reference functions are reported in Table |4j We chose Ac and Afitc to represent 
values computed form the level set approaches. The level sets no longer travel (on average) and 
hence have wave speeds with values close to zero. Note that the noise path is not the same as 
solving the SPDE for Table [3] and so we do not expect the values to be exactly the same, they 
are however close. The wave speed A^*^^ estimated by the minimization is identical for the two 
different reference functions (solving the SPDAE we do not have the same interpolation errors as 
when solving the SPDE). However, the choice of reference function is an issue for the SPDE and 
we consider this further in Section r3. 2. 2 1 

In Figure[2]we have plotted in (a) the space-time plot of solution of the SPDAE. The front starts 
at a; 200 and remains (on average) at that position throughout the computation illustrating 
that the wave does not travel (compare to Figure [l] (a)). In (b) for the two different initial data 
fco = l/v'2 (with mean 1.08522) and fco — 0.1 (with mean 2.74311) we have plotted the two 
distributions of the instantaneous wave speed A used to compute the wave speed through the 
minimization (with fc = l/y/2). Comparing with Figure [l] (b) we see similar distributions and 
greater variance with initial data with fco = 0.1 than fco — l/\/2 (as in Figure [I]). In (c) we 
have plotted for the two different initial data sets the instantaneous wave speed X{t) and the time 
averaged wave speed Aniin(i) with to = 0. We see faster convergence of the wave speed A(t) for 
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= 50, ti = 


100 


min 


7 


Ac 


Afitc 


ko 


= 1/V2,k 


= 1/V2 


1.08522 


1.08828 


2.421e-04 


-2.131C-04 


ko 


= l/V2,k 


= 0.1 


1.08522 


1.08828 


2.421e-04 


-2.131C-04 


ko 


= 0.1, fc = 


0.1 


2.74311 


2.76005 


-5.703e-02 


-6.404e-03 


ko 


= 0.1, fc = 


1/V2 


2.74311 


2.76005 


-5.703e-02 


-6.404e-03 



Table 4: Wave speeds computed from solving a single realization of the SPDAE with noise intensity 
= 0.1 and correlation length ^ = 0.1. 





A ■ or A^" 


A^ or A;(''^ 


Ac Afitc 


SPDE 
SPDAE 


1.08588 ± 0.19680 
1.08951 ± 0.19512 


1.08388 ± 2.73e-03 
1.08790 ± 2.39e-03 


1.08381 ±2.81e-03 1.08390 ±2.68e-03 
-4.0e-05 ± 2.0e-5 3.0e-05 ± 2.0e-5 



Table 5: Expected values of the wave speeds taken over 100 realizations solving the SPDE and 
the SPDAE. Initial data taken with ko ~ l/-\/2,and reference function with k = l/\/2. 



ko = 1/^2 than for ko = 0.1. 

(a) (b) (c) 




X 



Figure 2: (a) Space-time plot of a single realization of the frozen SPDE showing a travelling wave. 
In (b) distributions of the instantaneous wave speeds A computed for the two different initial data 
sets. In (c) we plot X{t) and the time averages Aniin(0 with to = 0, ti = t. 

Rather than looking at a single realization more physically meaningful results are found from 
taking the expectation over many realizations. In Table [5] we examine wave speeds based on 
100 realizations of both the SPDE and SPDAE. The different measures of the wave speed are in 
broad agreement. The larger uncertainty in Amin and A^^^j originates in the large variance in the 
instantaneous wave speeds A and is a drawback of the minimization approach. 

We also compare the profiles from the SPDE to profiles obtained from the SPDAE. To avoid 
the spreading of the wave we need to align individual realizations of the SPDE. We chose as a 
common reference the level set c(lOO). If we examine the final time profiles for the runs we find 
that the weak error, ||E {uspdae{WO)) — E (usp£)£;(100)) \\^2 for 10 realizations is « 0.0150 and 
for 100 realizations ~ 0.0144 and ~ 0.0117 with 1000 realizations. 

In Figure [s] we compare results for the SPDAE (a) and SPDE (b) for a range of different 
nonlinearity's a G {—1, —0.5, —0.3, 0, 0.3, 0.45} and noise intensities measured by fj,^ G [0, 1]. Initial 
data approximates a step function and the spatial correlation length of the noise ^ is that of the 
computational grid Ax. The results in (a), where Ax = C = 0-5 where agree with those in 
[21 12], reproduced in |16j . where the authors obtain a front velocity taking an average over an 
"appropriate time window" of Jj^u{x,t)dt and compare to a small noise analysis. In (b) we took 
a smaller spatial step Ax = 0.1 and have plotted the wave speed Amin computed both from 
minimization and from the level set, Ac on which the error bars are based. 
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We observe that the effect of the two different approximations to spatially white noise is to 
increase the speed of the wave. Note that some realizations where the wave is frozen in Figure |3] 



(a) fail to exist due to numerical instability, see Section 3.2.1 



2.5 



1.5 



(a) 







nr--1 - a 






' " " 










(b) 




Figure 3: Wave speeds Aj^in with increasing noise intensity for Stratonovich noise with correlation 
length equal to that of the grid. In (a) solving the SPDAE where the wave is frozen with ^ = 
0.5 = Ax (b) the SPDE with ^ = 0.1 = Ax. Each line corresponds to a different nonlinearity with 
a e {-1,-0.5,-0.3,0,0.3,0.45}. 



3.2.1 Numerical instability 

The numerical approximation of SDEs and SPDES where the solution is constrained in phase space 



is an area under development. For the Nagumo equations (14) (or (15)) u £ [0, 1], another typical 
example is a positivity constraint where u > 0. Numerical instability can lead to non-physical 
solutions and potentially to unphysical unbounded growth of the numerical solution. A number 
of approaches have been proposed to simulations to enforce constraints on the numerics and a 
review of these types of methods for SDEs is contained in [55]. One method to avoid unbounded 
growth in numerics from nonphysical solutions the nonlinearity and noise can be adapted as in 

[251 na El Eg. 



We found that solving the SPDEs (14) (or (15)) such instability was not an issue. However, 
when freezing the wave and solving the SPDAE ( |27[ ) did lead to non-physical solutions. In Figure|4] 
we have frozen the wave and show one realization at t = 34.7 (a) with an instantaneous wave speed 
10.98 and (b) t = 35.2 with an instantaneous speed —10.60. The non-physical regions where u < 
and u > 1 then grow in magnitude with further iterations. 

Figure [T] and Figure [2] show the distribution of the instantaneous wave speed A as computed for 



the SPDE and SPDAE respectively where the wave is frozen. The SPDAE system ( 27 ) includes the 
advection term Xvx , where A is a random variable with a particular distribution - which may lead 
to either large positive and/or negative values of A. Numerically this is particularly true for large 
noise intensities or small correlation lengths of the noise. The result of this is a loss of numerical 
stability. Although we were able to control unbounded growth by modifying the equation solved 
close to the u = and u = 1, direct comparisons of wave speeds to SPDE calculations showed this 
can lead to a bias in the estimate, so we do not include such results here. Hence, in Figure [3] (a) 
results for the SPDAE equation are reported with the expectation taken over solutions that existed 
to the final time. A large number of initial realizations was taken so that the final expectation is 
over at least 1000 realizations. Although we observe the same results calculating the wave speed 
based on level set methods, in general not taking the results where there is numerical blow up 
may bias the statistics. 
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Figure 4: Plot of realization of a single realization of the noise illustrating instability (a) a,tt = 34.7 
and (b) t = 35.2, with corresponding instantaneous wave speeds of 10.98 and —10.60. 



3.2.2 Choice of reference function it for the minimization 

We commented early that natural choice for the reference function would be to take either a 
deterministic travelling wave or the initial data. However, the width of the reference function it 
plays an important role in the computed wave speed for both the SPDE and SPDAE. If we take 
the reference function ii to be the Heaviside function then the minimization of the norm fails. 
We observe numerically that narrow reference functions can also lead to numerical failure of the 
minimization and indeed there may be more than one minimal position. 

To illustrate this we solve the SPDE with a = 0.25 and examine large noise intensity = 1 
combined with a small correlation length of f = 0.5. In Figure [s] (a) with u = Uk, k — l/\/2 
we see that the width of the computed front is larger than the width of the reference function u 
and in (b) is plotted the corresponding instantaneous wave speed X{t) - with time average 0.6483 
and variance of 6.4568. In (b) the reference function u = Uj^ has k = 0.1 and the width of the 
reference is larger than the solution. In (d) is plotted the corresponding instantaneous wave speed 
X{t) - with time average 0.6091 and variance over time of 3.8006. In (a) the computation of 
the minimization fails at a later time (T « 55) however in (b) computation was continued to 
T > 100. In (a) the minimization of the norm is dominated by the random fluctuations in the 
front which is avoided with a template function with larger support. Provided the width of the 
reference function is comparable or larger than the width of the front the computations are robust 
although convergence rates of the wave speed can be much slower for poor choices of the reference 
function. 

3.3 Effects of Stratonovich and Ito noise 

Accurate numerical calculations are notoriously difficult in the deterministic case when the wave 
profile depends on the leading profile of the wave, see for example |13] for the Nagumo equation or 
|29j for the Fisher equation. We consider from now on initial data that converges to the minimum 
speed wave in the deterministic case and take initial data uq = Uk close to a step function with 
k — 50. We examine the effects on wave speed and support of the front from changing the noise 
intensity and correlation length for both Stratonovich and Ito noise. 

First we examine the effects of Stratonovich noise on the travelling wave in the Nagumo 
equation. Figure |6] shows wave speed as noise intensity /i increases for four different corre- 
lation lengths £, = 0.1,0.5,1 and 10. On each plot are plotted different nonlinearities a — 
0.3, 0.25, 0, —0.25, —0.3, —0.5, —1. Each point on the plot is an average over 100 realizations and 
wave speeds measured both from minimization Amin and from the level set Ac. In Figure|6]we have 
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Figure 5; One realization of the solution with two different reference functions ii = Uk at time 
t = 50 (a) and (c) and the corresponding different distributions of the instantaneous wave speeds 
X{t) (b) and (d). In (a) and (b) k — 1/V2 and (c) and (d) k — 0.1. Note the smaller variance in 
(d) with k = 



plotted the corresponding average widths. We see that increasing the noise intensity increases the 
wave speed, where as increasing the correlation length of the noise decreases the wave speed. The 
two effects essentially cancel each other in (d) and we see no overall effect on the noise intensity on 
the wave speed. We can also examine the form of the wave profile. In Figure[6]we have plotted the 
corresponding average widths of the wave as noise intensity fi and correlation length ^ are changed. 
For large noise the width of the waves increase and this effect is again reduced as the correlation 
length is increased. For a spatial correlation length ^ = Aa; = 0.1 we have an approximation of 
white noise in space, for this case we see that for a = 0.45 and a — 0.25 the width of the wave 
increases and a larger computational domain is required. 

For Ito noise the effect of the noise on wave speed and width of the waves is less pronounced, 
see Figure |8] for the wave speed and Figure [9] for the corresponding width of the waves. We see 
that for large noise, in contrast to the Stratonovich case, a slight drop in the wave speed for a 
correlation length ^ < 10. As we change the noise intensity we see (for most nonlinearities) a drop 
in the width of the wave - and so the front is steeper on average and the effect is more pronounced 
for shorter correlation lengths. 



3.4 Computations using averaged quantities 

In general computing wave profiles using averaged quantities leads to the wave being 'polluted' by 
the spread of the individual waves (see [Hj). 

In (32 1 we propose using an expected value of the instantaneous wave speeds for the SPDAE. 
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Figure 6: Wave speeds Amin and for increasing Stratonovich noise intensity and diflFerent spatial 
correlation lengths (a) ^ = 0.1, (b) ^ = 0.5, (c) ^ = 1 and (d) ^ = 10. Increasing the noise intensity 
increases the expected wave speed where as increasing the correlation length decreases the expected 
wave speed. 



We fix a spatial correlation of ^ = 0.5. For the SPDAE if we solve with [i = 0.1, u ~ with 
k = 0.1 and v? — u^g, fco = l/v^ and 100 realizations then we obtain an estimate of a wave speed 



of 1.086. This compares with A„iin = 1.086 and A^ = 1.084 from solving the SPDE (14). If we 



examine the computed mean solution front we do not observe spreading of the wave front (see 



Figure 10 (a)). We also note from (b) that the distribution of X{t) has smaller variance than that 



from solving (27). For the SPDE we can implement a version (|32j) where we move the reference 
function using the expected values of the instantaneous wave speeds. In Figure [lO] (c) we plot the 
distribution of A for same parameters as in (b). The mean values agree although the distributions 
are different. In (d) we see that computed wave speeds using the average instantaneous speed and 
wave speeds A^ computed using the level set approach are the same over a range of nonlinearities. 
These are the same as those computed using the SPDE, compare to Figure |6](b). 

3.5 Additive noise 

We briefly consider the case of additive noise in the SPDE for which, unless the noise has some 
special properties, a solution will in general cease to exist at some finite time. 

We now change the parameter a in the nonlinearity to a = 0.1 and illustrate how the SPDAE 



approach deals with nucleation and extinction of waves. In Figure 12 we have plotted in (a) a single 



realization of the SPDE (so not frozen) showing nucleation and subsequent extinction {t « 98) of a 



Freezing Stochastic Travelling Waves 



16 



(a) (b) 




2 2 

mu mu 



Figure 7: Expected width of the wave for increasing Stratonovich noise intensity and different 
spatial correlation lengths (a) ^ = 0.1, (b) ^ = 0.5, (c) ^ = 1, and (d) ^ = 10. As the noise 
intensity is increased the expected width of the wave front increases where as for fixed intensity 
increasing the correlation length reduces the expected width. 



travelling wave. In (b) is plotted a single realization from computing using the SPDAE approach. 
We see the wave is fixed in the domain and at t w 50 a wave is nucleated at x « 100 by the additive 
noise. The computations are based on the original wave which remains fixed until it interacts with 
the nucleated wave and is annihilated at t « 94 when the computations stop when the wave cease 
to exist. In (c) and in (d) we have plotted mean profiles for the SPDE and the frozen SPDAE 
systems. In each case we see a well defined front from the averaging and individual nucleations 
and annihilations are no longer distinguishable (although in (d) a large solution pollutes the data 
at t w 130). 

4 Discussion 

We have examined level set based methods and minimization to a reference function methods to 
calculate the wave speed of a stochastic travelling wave. Our numerical results illustrate these give 
comparable results. Numerically we saw that for reference functions with support much smaller 
than the support of the travelling wave that the minimization may fail. Using the minimization 
technique for the SPDE (when it is not frozen) is more computationally expensive than the level 
set based methods as it requires interpolation at each time step. 

The algorithm described for freezing the wave and solving the SPDAE has several numerical 
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Figure 8: Expected wave speeds Amin and Ac for increasing Ito noise intensity and different spatial 
correlation lengths (a) ^ = 0.1, (b) ~ 0.5, (c) ^ = 1 and (d) — 10. As noise intensity is increased 
we see a slight drop in wave speed and little effect from the changing correlation length. 



advantages over simply solving the SPDE if the numerical instability issues could be over come. 
The frozen wave does not require a large computational domain for long time simulations and the 
generation of the noise path is not so computationally expensive. The cost of of the minimization 
when the wave is fixed is minimal as we simply need to compute two inner-products. However 
the advection term is nontrivial - and the loss of numerical stability is a real issue where some 
realizations fail to exist as ignoring results where there is numerical blow up may bias the statistics. 

Our investigation of the Nagumo equation has revealed interesting and new computational 
observations that we have not seen reported in the literature. Although it was known that for 
Stratonovich noise increasing noise intensity increases wave speed we have also seen it increases 
the support of the wave. In addition increasing the spatial correlation decreases the wave speed 
and decreases the support of the wave. The reverse is observed for Ito noise: the noise intensity 
seems to decrease the wave speed and correlation length has little influence on the speed decreases 
the support of the wave. 

For additive noise in the Nagumo equation we see that the wave speed is increased with the 
noise intensity like in the multiplicative case - this is probably because of the small perturbations 
ahead of the front that make the wave faster. 
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Figure 9: Expected width of waves for increasing Ito noise intensity and different spatial correlation 
lengths (a) ^ = 0.1, (b) ^ = 0.5, (c) ^ = 1 and (d) £, = 10. Increasing noise intensity narrows the 
width of the wave and this effect is mitigated by increasing the correlation length. Legend for all 
four plots is given in (d). 
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niose intensities. 
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Figure 12: Nucleation of travelling waves and annihilation for the Nagumo equation with a ~ 0.1. 
In (a) the space-time plot shows computations of the SPDE (no freezing) and in (b) the SPDAE 
where the wave is frozen. In (c) and (d) are plotted means over realizations for the frozen and 
travelling cases. 
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